v20200426.2342
Para executar os exemplos aqui apresentados, sugerimos que crie um projeto e baixe os seguintes arquivos para a pasta do mesmo:
Há diversas funções e procedimentos envolvidos neste texto. Como roteiro orientador, verifique:
Delineamento entre participantes
teste de independência / dissociação
Delineamento intraparticipantes
teste de concordância entre dois métodos de avaliação
\[~\]
|
1 Exposição e Desfecho |
\[~\]
|
2 Não avalia causa e efeito: apenas associação. |
“The public impression of scientific and medical uncertainty that resulted was a crucial element in maintaining the market of current smokers as well as recruiting new ones. Industry literature, for example, frequently pointed to the fact that nonsmokers could also develop lung cancer. Therefore, they argued, how could one attribute lung cancer to cigarette smoking?”
medical knowledge was always provisional and contingent. Just as drugs deemed “effective” do not work in every case, so too a cause of disease does not always result in disease. As Richard Doll would later explain, the epidemiologists had identified a cause of lung cancer (and other diseases), not the cause.
Another vocal critic of the lung cancer findings was Sir Ronald Fisher, the leading biometrician and geneticist in Great Britain during the first half of the twentieth century and a man deeply committed to bringing statistical analysis to genetics and agricultural experimentation. His 1925 book, Statistical Methods for Research Workers, quickly became a classic, leading to academic appointments at University College London and Cambridge University. Fisher’s critiques were similar to Berkson’s. The ethical impossibility of conducting a randomized experiment led him to question the results of the epidemiological studies. As a believer in genetic* notions of cancer causality, Fisher speculated that there was some constitutional factor that led individuals both to become smokers and to get lung cancer, even though smoking and lung cancer might not be causally related. Doll and Hill repeatedly rebutted this theory, returning to the critical question of how to account for the rise in lung cancers during the twentieth century if the disease was simply “constitutional.”
* 1923: One of the [A. Bradford] Hill’s innovations was the first randomized, double-blind clinical trial, designed to reduce investigator bias in the evaluation of clinical outcomes. […] This method, which drew on Fisher’s agricultural experimentation in genetics, became a critical new tool for evaluating medical treatments.
By the 1930s, the relationship of smoking to cancer was a topic of uresolved debate. It was the life insurance industry, which like the tobacco corporations grew by leaps and bounds in the first half of the twentieth century, that took the lead in understanding the effects of smoking on health.
By the early 1950s, however, it was abundantly clear that the evidence implicating cigarette smoking as a risk to health was now of a different order. First, the link between smoking and disease was categorical, outside the realm of individual clinical judgment.
They [tobacco industry] responded with a new and unprecedent public relations strategy. Its goals was to produce and sustain scientific skepticism and controversy in order to disrupt the emerging consensus on the harms of cigarette smoking. This strategy required intrusions into scientific process and procedure. […] So long as there appeared to be doubt, so long as the industry could assert “not proven”. […] The future of the cigarette would now depend on the successful production of a scientific controversy.
Quando alguém lhe traz uma opinião coincidente com sua crença, você a toma como um fato. Quando alguém lhe traz um fato discordante de sua crença, você o toma como uma opinião.
|
3 Tomada de decisão: Hipótese nula (\(H_0\)) e hipótese alternativa (\(H_1\)) |
|
Adotaremos, neste texto, a seguinte convenção:
|
Os testes baseados na estatística qui-quadrado podem ser divididos em dois tipos:
Teste de aderência uma variável nominal com duas ou mais categorias: teste de frequências hipotetizadas
Teste de independência entre duas variáveis nominais com duas ou mais categorias
Antes de discutirmos os testes, vamos explorar um pouco a distribuição que utilizam. Diferente da distribuição normal que têm domínio de \([- \infty, + \infty]\), estas são distribuições assimétricas, iniciadas em zero, que tendem assintoticamente a zero com o valor crescente de \(\chi^2\), \([0, + \infty]\).
O comando para calcular a probabilidade acumulada entre zero e cada valor de qui-quadrado é:
dchisq(x, df, ncp = 0, log = FALSE)
onde x é o valor de \(\chi^2\), df são os graus de liberdade (degrees of freedom) e ncp é o parâmetro de não centralidade (non-centrality parameter = 0, por default, correspondente à distribuição centrada de \(H_0\); voltaremos ao ncp adiante). O parâmetro log, desligado por default é para a possibilidade de devolver o logarítmo do valor (não o exploraremos neste texto).
Pode-se gerar curvas desta distribuição com o seguinte código, mostrando seu aspecto entre 2 e 6 graus de liberdade:
chi2.valor <- seq(from=0, to=40, by=0.01)
graus.liberdade <- 2:6
for (g in 1:length(graus.liberdade))
{
p <- dchisq(chi2.valor, df=graus.liberdade[g])
plot(chi2.valor, p, type="l",
main=paste("Distribuição de Qui^2, df=",
graus.liberdade[g],sep=""),
xlab="Qui-quadrado", ylab="Densidade")
} ou, caso prefira ver todas as curvas no mesmo gráfico, sugere-se:
source("friendlycolor.R")
cor <- 3
padrao <- 1
leg.txt <- c()
leg.cor <- c()
chi2.valor <- seq(from=0, to=40, by=0.01)
graus.liberdade <- 2:6
for (g in 1:length(graus.liberdade))
{
p <- dchisq(chi2.valor, df=graus.liberdade[g])
if (g==1)
{
plot(chi2.valor, p, type="l", lty=padrao,
col=friendlycolor(cor), lwd=2,
main="Distribuições de Qui^2",
xlab="Qui-quadrado", ylab="Densidade")
} else
{
lines(chi2.valor, p, lty=padrao,
col=friendlycolor(cor), lwd=2)
}
# guarda para a legenda
leg.txt <- c(leg.txt, paste("df = ",graus.liberdade[g],sep=""))
leg.cor <- c(leg.cor, friendlycolor(cor))
cor <- cor+6
padrao <- padrao+1
}
legend("topright",
leg.txt,
col=leg.cor,
lwd=2,
lty=1:5,
box.lwd=0, bg="transparent") Permite descobrir se um conjunto de frequências observadas difere de um outro conjunto de frequências hipotetizadas.
Cento e dez pessoas foram solicitadas a manifestar suas preferências com respeito a chocolates. Queremos verificar se alguma das marcas é preferida em detrimento de outras.
Se não são (\(H_0\)), devemos esperar aproximadamente o mesmo número de pessoas em cada categoria.
Com certeza, não haverá exatamente o mesmo número de participantes em cada categoria, mas os números devem ter uma distribuição com razoável proximidade.
O número de pessoas que escolheu apenas uma entre quatro marcas diferentes está em chocolate.xlsx.
O teste está implementado em TesteQuiQuadradoAderencia_Chocolate.R:
# TesteQuiQuadradoAderencia_Chocolate.R
library(readxl)
cat("Teste de aderencia a distribuicao multinomial\n")
df_chocolate <- read_excel("chocolate.xlsx")
prmatrix(df_chocolate, rowlab = rep("",nrow(df_chocolate)))
out <- chisq.test(df_chocolate$Pessoas, p=df_chocolate$probH0,
simulate.p.value = TRUE, B = 1e6)
df <- chisq.test(df_chocolate$Pessoas, p=df_chocolate$probH0)$parameter
X2 <- out$statistic
p <- out$p.value
n <- sum(chisq.test(df_chocolate$Pessoas, p=df_chocolate$probH0)$observed)
V <- sqrt((X2/n)/df)
cat("X^2 =", X2 , "gl =", df, "p exato =", p, "\n")
cat("\nHeuristica de significancia: X^2/gl > 2:\n", X2/df, "\n")
cat("\nResiduos estandardizados:\n", out$stdres, "\n")
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("V de Cramer = ", V,"\n",sep="")
cat("(Grau ", gV, "de dessemelhanca entre as distribuicoes observada e hipotetizada\n",sep="")
Observe o código:
Este código gera a seguinte saída:
Teste de aderencia a distribuicao multinomial
Chocolate Pessoas probH0
"A" "20" "0.25"
"B" "60" "0.25"
"C" "10" "0.25"
"D" "20" "0.25"
X^2 = 53.63636 gl = 3 p exato = 9.99999e-07
Heuristica de significancia: X^2/gl > 2:
17.87879
Residuos estandardizados:
-1.651446 7.156264 -3.853373 -1.651446
V de Cramer = 0.4031556
(Grau intermediariode dessemelhanca entre as distribuicoes observada e hipotetizada
Com \(p=9.999 \cdot 10^-5\), rejeita-se \(H_0\) e, portanto, temos evidência de que a escolha das pessoas pelas marcas de chocolate não é homogênea. Observando a tabela, vemos que o desequilíbrio ocorreu pela maior preferência pelo chocolate B.
| Os graus de liberdade em uma tabela com \(L\) linhas e \(C\) colunas são dados por \[df = (L-1) \cdot (C-1)\] Portanto, uma tabela \(2\text{ x }2\), tem somente um grau de liberdade, e uma tabela \(4\text{ x }3\) tem seis graus de liberdade. |
As probabilidades relacionadas com \(H_0\) não precisam ser sempre as mesmas. Por exemplo, em genética, é comum precisarmos avaliar como é a herança de um genótipo a partir dos fenótipos observados.
Duas linhagens puras de plantas foram cruzadas: uma com pétalas amarelas e outra com pétalas vermelhas. A primeira geração, \(F_1\), tem pétalas cor de laranja. Então \(F_1\) é cruzada entre si, gerando \(F_2\) com 320 plantas, das quais 182 têm pétalas cor de laranja, 61 amarelas e 77 vermelhas.Há duas possibilidades para explicar este resultado, e queremos saber qual é o caso para o fenótipo das pétalas das flores desta espécie:
Neste caso, os números esperados para \(F_2\) são, respectivamente, 160, 80 e 80, seguindo as proporções 2:1:1.
Neste caso, os números esperados para \(F_2\) são, respectivamente, 180 , 60 e 80 (60+20), seguindo as proporções 9:3:(3:1).
A qual das duas hipóteses a herança da cor das flores desta espécie adere?
Implementamos uma solução em eiras.quiaderencia.R e TesteQuiQuadradoAderencia_Flores.R.
Observe o código com nova função declarada:
# eiras.quiaderencia.R
# Teste Qui-quadrado de aderencia
eiras.quiaderencia <- function (df_dados, B=1e6)
{
valores <- as.vector(unlist(df_dados[,1]))
probabilidades <- as.vector(unlist(df_dados[,2]))
# teste robusto (bootstrapping)
robusto <- chisq.test(valores, p=probabilidades,
simulate.p.value = TRUE, B = B)
# usando o teste convencional para obter os graus de liberdade
tradicional <- chisq.test(valores, p=probabilidades)
df <- tradicional$parameter
X2 <- robusto$statistic
p <- robusto$p.value
n <- sum(robusto$observed)
V <- sqrt((X2/n)/df)
# resultados
cat("----------------------------------------------------\n")
cat("Teste de aderencia a distribuicao multinomial\n")
prmatrix(df_dados, rowlab = rep("",nrow(df_dados)))
cat("X^2 =", X2 , "gl =", df, "p exato =", p, "\n")
cat("\nHeuristica de significancia: X^2/gl > 2:\n", X2/df, "\n")
cat("\nResiduos estandardizados:\n", robusto$stdres, "\n")
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("V de Cramer = ", V,"\n",sep="")
cat("(Grau ", gV, " de dessemelhanca entre as distribuicoes observada e hipotetizada\n",sep="")
cat("----------------------------------------------------\n")
}
Criar uma função e armazená-la em um arquivo separado facilita seu uso (o nome do arquivo tem o mesmo nome da função só por conveniência; poderia ter qualquer nome). Esta função, eiras.quiaderencia(), recebe um data frame com duas colunas e o número de iterações para o bootstrapping (com \(10^6\) por default).
Adaptar o código para outros casos torna-se muito mais simples. No código anterior, seria preciso alterar todas as ocorrências de df_chocolate e os nomes das colunas Pessoas e probH0. Agora, só é necessário modificar TesteQuiQuadradoAderencia_Flores.R:
# TesteQuiQuadradoAderencia_Flores.R
library(readxl)
source ("eiras.quiaderencia.R")
# pega o arquivo e mostra seu conteúdo
df_flores <- read_excel("flores.xlsx")
cat("\nDados:\n")
prmatrix(df_flores, rowlab = rep("",nrow(df_flores)))
# chama os testes estatísticos
cat("\nHipotese de codominancia:\n")
eiras.quiaderencia(df_flores[,c(2,3)], B=1e6)
cat("\nHipotese de epistasis recessiva:\n")
eiras.quiaderencia(df_flores[,c(2,4)], B=1e6)
Para funcionar, carrega-se a função na memória (com source()), lemos o arquivo de interesse flores.xlsx e chamamos eiras.quiaderencia() duas vezes, indicando as colunas a serem testadas:
Este código produz a seguinte saída:
Dados:
CorPetalas Numero Codominancia Epistasis
"laranja" "182" "0.50" "0.5625"
"vermelha" " 77" "0.25" "0.2500"
"amarela" " 61" "0.25" "0.1875"
Hipotese de codominancia:
----------------------------------------------------
Teste de aderencia a distribuicao multinomial
Numero Codominancia
182 0.50
77 0.25
61 0.25
X^2 = 7.65 gl = 2 p exato = 0.02224998
Heuristica de significancia: X^2/gl > 2:
3.825
Residuos estandardizados:
2.459675 -0.3872983 -2.452889
V de Cramer = 0.1093303
(Grau pequeno de dessemelhanca entre as distribuicoes observada e hipotetizada
----------------------------------------------------
Hipotese de epistasis recessiva:
----------------------------------------------------
Teste de aderencia a distribuicao multinomial
Numero Epistasis
182 0.5625
77 0.2500
61 0.1875
X^2 = 0.1513889 gl = 2 p exato = 0.9324391
Heuristica de significancia: X^2/gl > 2:
0.07569444
Residuos estandardizados:
0.2253745 -0.3872983 0.143223
V de Cramer = 0.01538002
(Grau minimo de dessemelhanca entre as distribuicoes observada e hipotetizada
----------------------------------------------------
Neste exemplo buscamos não rejeitar \(H_0\). Para o habitual \(\alpha=0.05\) estes resultados sugerem que a herança das cores das pétalas desta espécie é explicada por epistasis recessiva.
O teste Qui-quadrado de Pearson permite que se descubra se existe um relacionamento ou efeito de interação entre duas variáveis qualitativas nominais com duas ou mais categorias.
Cada observação da unidade observacional deve ser alocada em apenas uma categoria de cada uma das duas variáveis nominais.
Considere um estudo que investiga a efetividade dos capacetes de segurança de bicicleta na prevenção de traumas cranianos. Os dados consistem de uma amostra aleatória de 793 indivíduos envolvidos em acidentes ciclísticos durante um certo período. Deseja-se testar, com \(\alpha=0.05\), se o uso do capacete tem funcionado como um fator de proteção.
Formulamos as hipóteses nula e alternativa:
\(H_0: \text{não existe associação entre uso do capacete e trauma craniano}\)
\(H_1: \text{existe associação entre uso do capacete e trauma craniano}\)
\(\alpha = 5\%\)
Os dados são:
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
print (tabela) Trauma + Trauma -
Capacete + 17 138
Capacete - 130 508
Sob a hipótese nula (\(H_0\)), i.e., se não existe associação entre o uso de capacete e o trauma craniano nos acidentes, os valores esperados são:
chi2 <- chisq.test(tabela)
print(round(chi2$expected,3)) Trauma + Trauma -
Capacete + 28.733 126.267
Capacete - 118.267 519.733
É possível transformar em proporções, obtendo:
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
tabela_prop <- tabela/sum(tabela)
print (round(tabela_prop,3)) Trauma + Trauma -
Capacete + 0.021 0.174
Capacete - 0.164 0.641
e, sob \(H_0\), os valores esperados são:
soma_col1 <- sum(tabela_prop[,1])
soma_col2 <- sum(tabela_prop[,2])
soma_lin1 <- sum(tabela_prop[1,])
soma_lin2 <- sum(tabela_prop[2,])
tabela_propesp <- tabela_prop
tabela_propesp[1,1] <- soma_lin1*soma_col1
tabela_propesp[1,2] <- soma_lin1*soma_col2
tabela_propesp[2,1] <- soma_lin2*soma_col1
tabela_propesp[2,2] <- soma_lin2*soma_col2
print(round(tabela_propesp,3)) Trauma + Trauma -
Capacete + 0.036 0.159
Capacete - 0.149 0.655
Note que o teste em R é feito diretamente com os valores observados (pois o tamanho da amostra importa) e seus resultados foram armazenados, acima, na variável chi2. O tipo desta variável é dado por:
str(chi2)List of 9
$ statistic: Named num 6.7
..- attr(*, "names")= chr "X-squared"
$ parameter: Named int 1
..- attr(*, "names")= chr "df"
$ p.value : num 0.00964
$ method : chr "Pearson's Chi-squared test with Yates' continuity correction"
$ data.name: chr "tabela"
$ observed : 'table' num [1:2, 1:2] 17 130 138 508
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
$ expected : num [1:2, 1:2] 28.7 118.3 126.3 519.7
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
$ residuals: 'table' num [1:2, 1:2] -2.189 1.079 1.044 -0.515
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
$ stdres : 'table' num [1:2, 1:2] -2.7 2.7 2.7 -2.7
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
- attr(*, "class")= chr "htest"
sapply(chi2, typeof) statistic parameter p.value method data.name observed
"double" "integer" "double" "character" "character" "double"
expected residuals stdres
"double" "double" "double"
É uma estrutura complexa, da qual muitos valores podem ser extraídos quando quisermos utilizá-los dentro de um RScript. A parte em chi2$expected, que já utilizamos, contém a tabela com os valores esperados sob \(H_0\).
O resultado principal do teste é exibido pela variável como um todo:
print (chi2)
Pearson's Chi-squared test with Yates' continuity correction
data: tabela
X-squared = 6.7001, df = 1, p-value = 0.009641
A estatística de teste tem o valor \(X^2=6.7001\), com 1 grau de liberdade (\(df=1\)), correspondendo ao valor-p de \(0.009641\) (que será discutido adiante). Com estes valores podemos decidir pela rejeição ou não-rejeição de \(H_0\).
No exemplo de uma tabela 2x2 há somente 1 grau de liberdade. O gráfico da distribuição de qui-quadrado (para melhor visibilidade, computamos \(0 \le \chi^2 \le 8\), pois o valor de \(\chi^2\) vai de zero a \(\infty\), e truncamos o eixo \(y\) em 1.2) é:
chi2.valor <- seq(from=0, to=8, by=0.01)
graus.liberdade <- 1
p <- dchisq(chi2.valor, df=graus.liberdade)
plot(chi2.valor, p, type="l", ylim=c(0, 1.2),
xlab="Qui-quadrado", ylab="Densidade"
)|
Conceitualmente, o teste Qui-quadrado é um teste de proporções, no qual a estatística de teste (i.e., o valor de \(X^2\)) é obtida pela diferença quadrática entre as frequências absolutas observadas e esperadas, dada por \[X^2 = \sum{ {(o - e)^2} \over{e} }\]. ou \[X^2 = {{(ad-bd)^2(a+b+c+d)} \over {(a+b)(c+d)(b+d)(a+c)}}\] Em R, experimente calcular a estatística na “força bruta”:
O valor de \(X^2\) aumenta com a discrepância crescente entre os valores observados e esperados. Note, porém, que o valor calculado pelo código acima é diferente do que aparece em:
Isto acontece por causa da correção para continuidade de Yates, que será comentada adiante. Na “força bruta”, a correção pode ser implementada com:
Obviamente, nenhum destes cálculos manuais é necessário. A função chisq.test() já os executa e armazena no objeto que denominamos, neste exemplo, chi2.
|
Nosso interesse é saber a partir de qual o valor de \(X^2\) a discrepância entre os valores observados e esperados é menos provável que o valor de \(\alpha\) escolhido. Isto é encontrar o valor que separa a área de \(5\%\) sob a curva de distribuição em sua cauda direita. Este valor é chamado de qui-quadrado crítico, \(\chi^2_c\).
| Note que o raciocínio é bicaudal, buscando associação positiva ou negativa, mas a operacionalização do teste só utiliza a cauda superior. Isto acontece porque \(chi^2=0\) indica igualdade numérica (e, portanto, também estatística) entre os valores esperados e observados. Consequentemente, procuramos \(\alpha=0.05\) somente na cauda direita. |
Em R, podemos encontrar \(\chi^2_c\) pelo complemento (\(1-\alpha=0.95\)) com
chi2.critico <- qchisq(p=0.95, df=1)
cat("qui-quadrado crítico (df=1, alfa=5%): ", chi2.critico, "\n", sep="")qui-quadrado crítico (df=1, alfa=5%): 3.841459
pois este é o valor de \(\chi^2\) que deixa \(95\%\) da área sob a curva à sua esquerda.
Graficamente, podemos hachurar em azul a área correspondente a \(\alpha\) e traçar uma linha pontilhada vertical em \(\chi^2_c\) com:
source("friendlycolor.R")
# vetor com valores de 0 a 8
chi2.valor <- seq(from=0, to=8, by=0.01)
graus.liberdade <- 1
# vetor com probabilidades correspondentes
p <- dchisq(chi2.valor, df=graus.liberdade)
# exibe o gráfico
plot(chi2.valor, p, type="l", ylim=c(0, 1.2),
xlab="Qui-quadrado", ylab="Densidade")
# vetor com os valores acima de qui critico
chi2.critico <- qchisq(p=0.95, df=1)
chi2maioralfa <- chi2.valor[chi2.valor>chi2.critico]
# probabilidades correspondentes
pmaioralfa <- dchisq(chi2maioralfa, df=graus.liberdade)
# hachura sob a curva
chi2maioralfa <- c(min(chi2maioralfa), chi2maioralfa, max(chi2maioralfa))
pmaioralfa <- c(0, pmaioralfa, 0)
polygon(chi2maioralfa, pmaioralfa, col=friendlycolor(10), border = NA)
# reforca a linha
lines(chi2maioralfa[2:(length(chi2maioralfa)-1)],
pmaioralfa[2:(length(pmaioralfa)-1)],
col=friendlycolor(10), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2.critico, lty=2) Esta área hachurada é a região de rejeição de \(H_0\).
A estatística \(\chi^2\) é maior quanto maior for a discrepância entre os valores observados e esperados. Para um grau de liberdade e \(\alpha=0.05\), computamos \(\chi^2_c=\) 3.8414588. No exemplo de capacete e trauma, acima, rejeitamos \(H_0\) porque calculamos \(X^2=\) 6.7001206, valor maior que o valor crítico.
À direita do valor calculado \(X^2=\) 6.7001206 define-se uma área, que corresponde ao valor-p. Este valor também foi calculado quando aplicamos o teste e guardamos seu resultado em chi2, acima; o valor-p encontra-se acessível em chi2$p.value=0.0096406. Como \(p < \alpha\), concluimos pela rejeição de \(H_0\).
Podemos hachurar em laranja a área correspondente ao valor-\(p\) e traçar uma linha pontilhada vertical na altura do \(X^2\) calculado, adicione mais linhas ao código anterior, mas precisamos da variável chi2$statistic que retorna de chisq.tes(). O gráfico, agora, tem muitos elementos e, portanto, também devemos adicionar uma legenda.
O código completo é:
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
chi2 <- chisq.test(tabela)source("friendlycolor.R")
# vetor com valores de 0 a 8
chi2.valor <- seq(from=0, to=8, by=0.01)
graus.liberdade <- 1
# vetor com probabilidades correspondentes
p <- dchisq(chi2.valor, df=graus.liberdade)
# exibe o gráfico
plot(chi2.valor, p, type="l", ylim=c(0, 1.2),
xlab="Qui-quadrado", ylab="Densidade")
# vetor com os valores acima de qui critico
chi2.critico <- qchisq(p=0.95, df=1)
chi2maioralfa <- chi2.valor[chi2.valor>chi2.critico]
# probabilidades correspondentes
pmaioralfa <- dchisq(chi2maioralfa, df=graus.liberdade)
# hachura sob a curva
chi2maioralfa <- c(min(chi2maioralfa), chi2maioralfa, max(chi2maioralfa))
pmaioralfa <- c(0, pmaioralfa, 0)
polygon(chi2maioralfa, pmaioralfa, col=friendlycolor(10), border = NA)
# reforca a linha
lines(chi2maioralfa[2:(length(chi2maioralfa)-1)],
pmaioralfa[2:(length(pmaioralfa)-1)],
col=friendlycolor(10), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2.critico, lty=2)
# ########## Area a direita de p ############
# vetor com os valores acima de qui calculado
chi2maiorp <- chi2.valor[chi2.valor>chi2$statistic]
# probabilidades correspondentes
pmaiorp <- dchisq(chi2maiorp, df=graus.liberdade)
# hachura sob a curva
chi2maiorp <- c(min(chi2maiorp), chi2maiorp, max(chi2maiorp))
pmaiorp <- c(0, pmaiorp, 0)
polygon(chi2maiorp, pmaiorp, col=friendlycolor(20), border = NA)
# reforca a linha
lines(chi2maiorp[2:(length(chi2maiorp)-1)],
pmaiorp[2:(length(pmaiorp)-1)],
col=friendlycolor(20), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2$statistic, lty=3)
# legenda
legend("topright",
c("Qui^2 crítico","Qui^2 obs.", "alfa", "p"),
col=c("black", "black", friendlycolor(10), friendlycolor(20)),
lty=c(2, 3, 1, 1),
lwd=c(1, 1, 5, 5),
box.lwd=0, bg="white",
cex=0.8) |
As duas formas de decisão são igualmente válidas e equivalentes. Exibir o valor-p é a forma mais moderna e recomendada atualmente. |
Repare nas saídas da função chisq.test(), cujo cabeçalho indica:
Pearson’s Chi-squared test with Yates’ continuity correction
Esta é uma versão tradicional do teste, para a qual a correção de Yates é aplicável por default às tabelas 2x2 (mais detalhes serão vistos adiante). Podemos evitá-la com:
chi2ny <- chisq.test(tabela, correct=FALSE)
print(chi2ny)
Pearson's Chi-squared test
data: tabela
X-squared = 7.3099, df = 1, p-value = 0.006858
Além disto, há várias condições para a aplicação do teste Qui-quadrado (comentadas adiante) que podem, muitas vezes, restringir seu uso.
Uma alternativa é computar o teste em sua versão robusta. A função chisq.test() consegue simular a distribuição por bootstrapping, alterando-se o comando para:
chi2.robusto <- chisq.test(tabela, correct=FALSE,
simulate.p.value = TRUE, B=1e6)
print(chi2.robusto)
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: tabela
X-squared = 7.3099, df = NA, p-value = 0.007638
Esta versão é simulada com \(1e6 = 1 \cdot 10^6 = 1000000\) iterações.
| Como a amostra é grande, os valores de \(X^2\) e \(p\) não mudam muito entre as diferentes versões. Há um pequeno bug na implementação da versão robusta, que não calcula os graus de liberdade (sabemos que \(df=1\)). |
A decisão, neste caso, não mudou, tanto pelo critério do \(\chi^2_c\) quanto pelo \(p\): rejeita-se \(H_0\).
|
Odds e probabilidade Odds (traduzido habitualmente por chance) e probabilidade são formas equivalentes para expressar possibilidades e relacionadas por: \(\textit{Odds} = {{\text{Probabilidade}} \over{1-\text{Probabilidade}}}\) e \(\text{Probabilidade} = {{\textit{Odds}} \over{1+\textit{Odds}}}\) Por exemplo, uma probabilidade de \(80\%\) corresponde a \(\textit{Odds} = {{0.8} \over {1-0.8}} = {0.8 \over 0.2} = 4\) indicando que dizer “\(80\%\) de probabilidade” é equivalente a dizer “4 vezes mais chance de ocorrer do que não ocorrer” um determinado evento. Resersamente, Odds de 2 é: \(\text{Probabilidade} = {2 \over {1+2}} = {2 \over 3} \approx 66.67\%\) e \(\textit{Odds}=1\) resulta em: \(\text{Probabilidade} = {1 \over {1+1}} = {1 \over 2} = 50\%\) Associa-se o máximo de incerteza com probabilidade de \(50\%\); por este motivo, \(Odds=1\) será um valor necessário para decisões estatísticas, adiante. |
Há dois tipos principais de risco relativo: razão de riscos (RR, risk ratio) e razão de chances (OR, odds ratio).
Vamos considerar a seguinte tabela:
Tem efeito Não tem efeito
Tem exposição a b a+b
Não tem exposição c d c+d
a+c b+d total
Esta tabela de contingência genericamente relaciona exposição (e.g., hábito de fumar) com o efeito (e.g., câncer de pulmão). Traz as contagens dos:
Riscos são probabilidades. A razão de riscos é dada por:
\(RR = \LARGE{{{{a} \over {a+b}} \over {{{c} \over {c+d}}}}}\)
Repare que \({a} \over {a+b}\) é a probabilidade de um indivíduo exposto ter o efeito (e.g., fumante ter câncer); \({c} \over {c+d}\) é a probabilidade de um indivíduo não exposto ter o efeito (e.g., não-fumante ter câncer).
RR, portanto, é um odds: quantas vezes mais é esperado o efeito (ocorrência de câncer) de quem é exposto (fumante) em comparação com quem não é exposto (não fumante);
Odds ratio como diz o nome, é uma razão entre dois odds, dada por:
\(OR = \LARGE{{ { { {a} \over {a+c} } \over { {c} \over {a+c} } } \over { { {b} \over {b+d} } \over { {d} \over {b+d} } } }}\)
Considere, agora, que:
Interessante, também, é que (com alguma álgebra) \(OR\) reduz-se a:
\(OR = \LARGE{ {{a \cdot d} \over {{b \cdot c}} } }\)
fácil de lembrar e assim conhecido como “produto cruzado” (com o defeito de ocultar o motivo pelo qual é um odds ratio). Por outro lado, como produto cruzado, percebe-se que é uma medida de quanto pesa a:
em relação à
| Os valores de \(RR\) e \(OR\) são próximos quando menos que \(5\%\) dos não expostos têm efeito. |
Em R, o cálculo de \(OR\) aparece no Teste Exato de Fisher. No exemplo de capacete e trauma, temos:
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
print(tabela)
ft <- fisher.test(tabela) # Teste de OR robusto
print(ft) Trauma + Trauma -
Capacete + 17 138
Capacete - 130 508
Fisher's Exact Test for Count Data
data: tabela
p-value = 0.005688
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.2630993 0.8357342
sample estimates:
odds ratio
0.4817645
A estimativa pontual de odds ratio está em ft$estimate=0.4817645, mas para a decisão é obrigatório verificar seu intervalo de confiança 95% (\(IC95\)), com limites inferior e superior respectivamente guardados em ft$conf.int[1]=0.2630993 e ft$conf.int[2]=0.8357342.
Como a ausência de efeito é dada por \(OR=1\) (correspondendo a \(H_0\)) e este valor unitário está fora do intervalo, rejeita-se \(H_0\). Além disto, o \(IC95\) é menor que o esperado por \(H_0\) e, portanto, o uso de capacete está associado com redução dos traumas cranianos.
|
Caso a tabela fosse construída com as colunas em posições trocadas, obter-se-ia o seguinte:
O odds ratio é agora maior que o valor unitário, o qual está fora do intervalo, rejeitando-se a hipótese nula e a associação é positiva: a conclusão é a mesma, o uso de capacete está associado com o aumento de “Trauma -”, portanto é protetor. Aliás, observe que
são os mesmos valores obtidos anteriormente.
|
Há outros pacotes em R que calculam risk ratio e odds ratio, como por exemplo:
library(epitools)
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
epitools::oddsratio(tabela, method="fisher")que produz os mesmos valores para a estimativa pontual e intervalo de confiança que o teste exato de Fisher:
$data
Trauma + Trauma - Total
Capacete + 17 138 155
Capacete - 130 508 638
Total 147 646 793
$measure
NA
odds ratio with 95% C.I. estimate lower upper
Capacete + 1.0000000 NA NA
Capacete - 0.4817645 0.2630993 0.8357342
$p.value
NA
two-sided midp.exact fisher.exact chi.square
Capacete + NA NA NA
Capacete - 0.005088331 0.005688151 0.006857643
$correction
[1] FALSE
attr(,"method")
[1] "Conditional MLE & exact CI from 'fisher.test'"
Há mais três métodos disponíveis, de acordo com a documentação da função.
A library epitools também permite o cálculo de \(RR\) com intervalo de confiança. Ainda que não seja adequada para o presente exemplo, verifiquem sua sintaxe. Há três métodos disponíveis, como por exemplo:
library(epitools)
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
epitools::riskratio(tabela, method="wald")$data
Trauma + Trauma - Total
Capacete + 17 138 155
Capacete - 130 508 638
Total 147 646 793
$measure
NA
risk ratio with 95% C.I. estimate lower upper
Capacete + 1.0000000 NA NA
Capacete - 0.8943256 0.8357184 0.9570428
$p.value
NA
two-sided midp.exact fisher.exact chi.square
Capacete + NA NA NA
Capacete - 0.005088331 0.005688151 0.006857643
$correction
[1] FALSE
attr(,"method")
[1] "Unconditional MLE & normal approximation (Wald) CI"
e também uma variante com bootstrapping:
library(epitools)
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
epitools::riskratio(tabela, replicates=1e6)$data
Trauma + Trauma - Total
Capacete + 17 138 155
Capacete - 130 508 638
Total 147 646 793
$measure
NA
risk ratio with 95% C.I. estimate lower upper
Capacete + 1.0000000 NA NA
Capacete - 0.8943256 0.8357184 0.9570428
$p.value
NA
two-sided midp.exact fisher.exact chi.square
Capacete + NA NA NA
Capacete - 0.005088331 0.005688151 0.006857643
$correction
[1] FALSE
attr(,"method")
[1] "Unconditional MLE & normal approximation (Wald) CI"
| Lembre-se de nunca decidir a respeito de risk ratio ou odds ratio sem considerar o intervalo de confiança. Por isso, evite fazer o cálculo manualmente: use R! |
|
4 Simulação e Gráficos |
Em tabelacontingencia2x2_capacetetrauma.R, que você pode adaptar para as suas necessidades, processamos o mesmo exemplo (procure a função sink() que desvia o resultado da tela para um arquivo texto). Além da estatística \(X^2\), computamos outros valores interessantes. Observe sua saída (tabelacontingencia2x2_capacetetrauma.txt):
source("tabelacontingencia2x2_capacetetrauma.R")
Resultados armazenados em tabelacontingencia2x2_capacetetrauma.txt
Em Qui2.R também processamos o mesmo exemplo. Esta versão aplica bootstrapping de um jeito mais manual e produz gráficos associados aos conceitos que discutimos.
Execute com:
source("Qui2.R")
**** TABELA DE CONTINGÊNCIA ****
Considere:
Coluna+ Coluna-
Linha + a b Total lin. +
Linha - c d Total lin. -
Total col. + Total col. - Total geral
Informe:
Nome nas linhas (ate 10 letras): Capacete
Nome nas colunas (ate 10 letras): Trauma
Definido:
Trauma + Trauma -
Capacete + a b Total de Trauma +
Capacete - c d Total de Trauma -
Total de Capacete + Total de Capacete - Total geral
Forneca os valores de a, b, c, d:
a: 17
b: 138
c: 130
d: 508
Defina alfa = prob. erro do tipo I).
(número entre 0 e 1, default = 0.05).
alfa: 0.05
Quantas iteracoes para simular?
(entre um numero inteiro; default n=1e4)
iteracoes: 1e4
Exibir tabelas? (exibir lentifica a simulacao)
0=nao, 1=sim; default eh 0: 0
Observe a simulação em andamento na área de Plots. Quando terminar, procure os gráficos finalizados. Os principais são:
| O valor de \(\beta\) e o poder associados, aqui, são os valores a posteriori, portanto inúteis para qualquer decisão |
|
Flashback: Recorde a tomada de decisão com a distribuição binomial: |
|
Odds ratio é uma medida da intensidade de associação entre duas variáveis nominais dicotômicas de exposição e desfecho. De acordo com SHARPE, D. (2015) Your chi-square test is statistically significant: now what? Practical Assessment, Research & Evaluation, 20(8). É, também, uma medida de tamanho de efeito (não depende do tamanho da amostra, é adimensional e adirecional) e sua intensidade pode ser classificada em:
|
|
O programa já indica o tamanho de efeito de acordo com
Tamanho de efeito:
|
O teste Qui-quadrado opera em tabelas maiores, com \(L\) linhas e \(C\) colunas.
Por exemplo, três quimioterápicos foram comparados quanto ao efeito colateral (náusea). Considere \(\alpha=0.05\). Há diferença entre as drogas?
Os dados obtidos foram:
Nausea Sem_nausea
Droga A 3 5
Droga B 7 2
Droga C 6 3
Em Qui2_LxC.R implementamos este exemplo (versão robusta):
source("Qui2_LxC.R") ComNausea SemNausea
A 3 5
B 7 2
C 6 3
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: tcrxc
X-squared = 3.0559, df = NA, p-value = 0.2744
V de Cramer = 0.3428334
Observe que não rejeita-se \(H_0\). Os dois critérios são sempre equivalentes:
Abra o código para ver como foi implementado este teste.
|
5 Condições para o teste Qui-quadrado (métodos não robustos) |
Antes das implementações computacionais e a facilidade do uso do R, a partir de uma tabela 2x2 era necessário construir a tabela dos valores esperados.
Trauma + Trauma -
Capacete + 17 138 155
Capacete - 130 508 638
147 646 793
Trauma + Trauma -
Capacete + 28.7 126.3
Capacete - 118.3 519.7
Cada valor esperado é obtido pelo produto das marginais, dividido pelo total geral:
Avaliando os dados observados e os valores esperados, era necessário verificar se o teste Qui-quadrado podia ser aplicado. Eram as seguintes condições:
SIEGEL, S. & CASTELLAN Jr., N.J. (1988) Nonparametric statistics for the behavioral sciences. 2ª ed. NY: McGraw-Hill, p. 123:
COCHAN, W.G. (1954) Some Methods for Strengthening the Common \(\chi^2\) Tests. Biometrics: 10(4): 417-451:
Então, o valor de \(X^2\) era dado por
\[X^2 = \sum{ ({\text{observado}-\text{esperado}})^2 \over {\text{observado}} }\]
neste exemplo calculado por:
\(X^2 = { {{(17-28.7)^2} \over {28.7}} + {{(138-126.3)^2} \over {126.3}} + {{(130-118.3)^2} \over {118.3}} + {{(508-519.7)^2} \over {519.7}} } \approx 7.31\)
ou, em tabelas 2x2, era recomendado usar a correção de continuidade de Yates, dada por:
\[X^2 = \sum{ (|{\text{observado}-\text{esperado}}|-0.5)^2 \over {\text{observado}} }\]
Neste exemplo:
\(X^2 = { {{(|17-28.7|-0.5)^2} \over {28.7}} + {{(|138-126.3|-0.5)^2} \over {126.3}} + {{(|130-118.3|-0.5)^2} \over {118.3}} + {{(|508-519.7|-0.5)^2} \over {519.7}} } \approx 6.7\)
Então, era necessário consultar uma tabela com os valores críticos de \(\chi^2_c\) previamente calculados para a tomada de decisão (sem obter, portanto, o valor-\(p\)). A tabela permitia escolher o valor de \(\alpha\) e, sabendo-se os graus de liberdade (neste caso \(\nu = 1\)), localizar \(\chi^2_c\) (neste exmeplo igual a \(3.841\)):
Caso \(\chi^2\) não atendesse as exigências para ser aplicado, a alternativa (se fosse tabela 2x2) era o Teste Exato de Fisher. Este envolve calcular o valor-\(p\) por fatoriais de tabelas progressivamente mais extremas. Por exemplo:
Neste caso o valor exato de \(p\) era calculado e comparado diretamente com \(\alpha\) para a tomada de decisão.
Esta análise permite, quando se rejeita \(H_0\), localizar quais células da tabela de contingência mais contribuíram para esta rejeição.
Para descobrir se existe uma relação entre tabagismo e etilismo a partir de um estudo realizado com 100 estudantes universitários, encontrou-se:
| Tabagista | Não tabagista | |
|---|---|---|
| Etilista | 50 | 15 |
| Não etilista | 20 | 25 |
|
A manipulação de dados em R é poderosa, mas nem sempre fácil. Para a entrada de dados na forma de uma matriz 2x2 o seguinte código funciona: Obtendo-se:
O problema está na entrada de dados. Somente espaços em branco podem ser usados para a montagem da tabela. O alinhamento é visual, mas não necessário para o R. Funciona igualmente
Portanto, pode ser ruim a visualização de sua entrada. Uma opção mais “limpa” é sempre colocar os dados em planilhas no formato .xls ou .xlsx e usar a função readxl::read_excel() para utilizá-los em R:
A função read_excel() avisa que teve que suprir um nome (por causa da célula A1 que está vazia), mas esta célula será desprezada adiante. No entanto, temos um problema a resolver. A função read_excel() produz um data frame e a função chisq.test() precisa de uma matriz para sua entrada:
Temos, portanto, que converter TCnew para a representação em matrix. Existe uma função para isto:
Ainda não é o desejado. TCmatrix tem 2 linhas mas 3 colunas, além de ter perdido os nomes das linhas que TCnew trazia em sua primeira coluna:
A solução, portanto, é transferir o conteúdo de TCnew como nomes das linhas de TCmatrix e, então, deletar a coluna 1 de TCmatrix que contém valores NA.
Colocando tudo junto, para terminar com uma matriz chamada TC (como era originalmente), podemos ler a planilha em uma variável (e.g., TCtmp), reorganizar seu conteúdo em TC e remover TCtmp da memória com a função remove(). O código enxuto é:
Uma mensagem:
apareceu aqui. Warnings são avisos do R que não impedem o funcionamento do programa (avalie, caso a caso, se precisa corrigir algo; aqui não é necessário). Adiante mostramos como é removê-la. |
Neste exemplo, o objetivo é executar o teste Qui-quadrado mas, no caso de rejeitar-se \(H_0\), executar uma análise post hoc. Observe o código de TabelaContingencia2x2_TabagismoEtilismo.R que pega os dados da planilha tabagismo_e_etilismo_2x2.xlsx:
# TabelaContingencia2x2_TabagismoEtilismo.R
# desabilita warnings
options(warn=-1)
# Tabela de Dancey&Reidy (2019, p. 275)
TCtmp <- read_excel("tabagismo_e_etilismo_2x2.xlsx")
TC <- data.matrix(TCtmp)
rownames(TC) <-as.character(unlist(TCtmp[1:2,1]))
TC <- TC[,-1]
remove(TCtmp)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
print(TC)
cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05
cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(TC,simulate.p.value=TRUE,B=1e6))
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=1), "\n", sep="")
N <- sum(res$observed)
nL <- nrow(TC) # ou dim(TC)[1]
nC <- ncol(TC) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
cat("\nResiduos ajustados standardizados corrigidos por momento (MCSTARs):\n")
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
print(MCSTAR)
alfaBonf <- alfa/df
zcrit <- abs(qnorm(alfaBonf/2))
cat("\n|MCSTAR critico| (alfaBonferroni=5%/",df,") = ",zcrit,"\n", sep="")
cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
V <- sqrt(phi2/Dim) # V ou C de Cramer
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("\nV de Cramer =", V, "\nGrau", gV,
"de dependencia entre as duas variaveis nominais\n")
cat("\nTeste de razao de chances (OR) robusto:\n")
resft <- fisher.test(TC,conf.level = 1-alfa) # Teste de OR robusto
print(resft)
# reabilita warnings
options(warn=0)
que produz a seguinte saída:
New names:
* `` -> ...1
------------------------------------------
Dados:
------------------------------------------
Tabagista NaoTabagista
Etilista 50 15
NaoEtilista 20 25
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: TC
X-squared = 12.121, df = NA, p-value = 0.000666
X^2 critico de 95% = 3.841459
Graus de liberdade (nao fornecidos por bootstrapping): 1
Heuristica de significancia (rej. H0 se X^2/gl > 2): 12.12149
------------------------------------------
Analise post hoc:
------------------------------------------
Residuos ajustados standardizados corrigidos por momento (MCSTARs):
Tabagista NaoTabagista
Etilista 6.963186 -6.963186
NaoEtilista -6.963186 6.963186
|MCSTAR critico| (alfaBonferroni=5%/1) = 1.959964
------------------------------------------
Analise de significancia pratica:
------------------------------------------
V de Cramer = 0.3319569
Grau intermediario de dependencia entre as duas variaveis nominais
Teste de razao de chances (OR) robusto:
Fisher's Exact Test for Count Data
data: TC
p-value = 0.0006368
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.694576 10.331999
sample estimates:
odds ratio
4.107342
A novidade está na sessão Residuos ajustados standardizados corrigidos por momento (MCSTARs). Analise pelos valores positivos acima do valor assinalado como [MCSTAR critico]: são escores \(z\) relativos a quanto cada célula contribuiu para o qui-quadrado estimado. Neste exemplo, a diagonal principal é a responsável pela rejeição de \(H_0\), indicando associação concordante entre etilismo e tabagismo.
|
A única “sujeira” deste código é:
Isto desaparece se colocarmos algo, que será desprezado, na célula A1 da planilha. Experimente rodar o código com a planilha tabagismo_e_etilismo_2x2_v2.xlsx. Note, também, o uso das linhas:
|
Segundo este autor, a correção de Yates não deve ser usada:
Note a relação entre esta forma de calcular o qui-quadrado e o odds-ratio:
O teste qui-quadrado tradicional trata os dois fatores com categorias nominais. Há muitos casos, porém, em que estes fatores são ordinais. Embora não seja incorreto utilizar o qui-quadrado (toda variável ordinal é também nominal), alguma informação está sendo perdida.
Existe uma medida, o Gama de Goodman-Kruskal, que considera ambas as variáveis como categóricas ordinais, potencialmente melhorando o resultado do teste. É uma estatistica nao-parametrica que mede a força de associação em dados cruzados em tabela, com as hipóteses:
\[H_o: \gamma = 0\] \[H_o: \gamma \ne 0\]
Neste estudo, Goodman (1979) desenvolve a análise com variáveis ordinais.
Os dados desta tabela são, originalmente, de:
Duncan, O.D. (1979), How Destination Depends on Origin in the Occupational Mobility Table. American Journal of Sociology, 84, 793-803.
As categorias são as seguintes:
Portanto, a variável de status ocupacional é, claramente, qualitativa ordinal.
| O repositório do R tem mais que pacotes. Há vários bancos de dados armazenados nele, para uso em exemplos. Aqui usaremos os dados deste estudo, disponíveis no pacote “dataset”. |
Implementado em TCrXc_StatusOcupacional.R:
# TCrXc_StatusOcupacional.R
library(datasets)
library(DescTools)
eiras.show.MCSTAR <- function (MCSTAR, zcrit)
{
MCSTARtxt <- round(MCSTAR,3)
for (r in 1:nrow(MCSTAR))
{
for (c in 1:nrow(MCSTAR))
{
if (MCSTAR[r,c] > zcrit)
{
MCSTARtxt[r,c] <- paste("#",MCSTARtxt[r,c],sep="")
} else
{
MCSTARtxt[r,c] <- paste(" ",MCSTARtxt[r,c],sep="")
}
}
}
return(MCSTARtxt)
}
# desabilita warnings
options(warn=-1)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
# Duncan, O.D. (1979), How Destination Depends on Origin in the
# Occupational Mobility Table. American Journal of Sociology 84, 793-803.
tcrxc <- datasets::occupationalStatus
# para recolocar as categorias originais
categories <- c("I","II","III","IV","Va","Vb","VI","VII")
rownames(tcrxc) <- categories
colnames(tcrxc) <- categories
prmatrix(tcrxc)
cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05
cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(tcrxc,simulate.p.value=TRUE,B=1e6))
N <- sum(res$observed)
nL <- nrow(tcrxc) # ou dim(TC)[1]
nC <- ncol(tcrxc) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=df), "\n", sep="")
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
cat("\nResiduos ajustados standarizados corrigidos por momento (MCSTARs):\n")
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
alfaBonf <- alfa/df
zcrit <- abs(qnorm(alfaBonf/2))
print(eiras.show.MCSTAR(MCSTAR,zcrit))
cat("\n|MCSTAR critico| (alfaBonferroni=5%/",df,") = ",zcrit,"\n", sep="")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
V <- sqrt(phi2/Dim) # V ou C de Cramer
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("\nV de Cramer =", V, "\nGrau", gV,
"de dependencia entre as duas variaveis nominais\n")
# Goodman, L. A. (1979) Simple Models for the Analysis of Association
# in Cross-Classifications having Ordered Categories.
# J. Am. Stat. Assoc., 74 (367), 537–552.
gama <- DescTools::GoodmanKruskalGamma(tcrxc)
cat("\nGama de Goodman-Kruskal = ", gama,"\n", sep="")
# EDWARDES, MD & BALTZAN, M (2000) The generalization of the odds ratio,
# rik ratio and risk difference to rxk tables.
# Statistics in Medicine 19:1901-14.
ORg <- (1+gama)/(1-gama)
cat("\nOR generalizado = (1+",gama,")/(1-",gama,") = ", ORg,"\n", sep="")
# reabilita warnings
options(warn=0)
A saída é:
------------------------------------------
Dados:
------------------------------------------
I II III IV Va Vb VI VII
I 50 19 26 8 7 11 6 2
II 16 40 34 18 11 20 8 3
III 12 35 65 66 35 88 23 21
IV 11 20 58 110 40 183 64 32
Va 2 8 12 23 25 46 28 12
Vb 12 28 102 162 90 554 230 177
VI 0 6 19 40 21 158 143 71
VII 0 3 14 32 15 126 91 106
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: tcrxc
X-squared = 1416, df = NA, p-value = 1e-06
X^2 critico de 95% = 66.33865
Graus de liberdade (nao fornecidos por bootstrapping): 49
Heuristica de significancia (rej. H0 se X^2/gl > 2): 28.89877
------------------------------------------
Analise post hoc:
------------------------------------------
Residuos ajustados standarizados corrigidos por momento (MCSTARs):
destination
origin I II III IV Va Vb VI VII
I #28.022 #6.466 #4.851 -2.711 -0.804 -7.091 -4.336 -4.284
II #6.535 #15.194 #6.477 -0.475 0.201 -6.217 -4.43 -4.437
III 0.706 #6.01 #7.195 #3.979 2.782 -3.966 -6.129 -4.134
IV -1.369 -0.926 1.7 #6.772 0.826 0.847 -3.453 -5.132
Va -1.436 0.409 -0.87 0.701 #5.188 -1.363 0.388 -1.982
Vb -6.546 -6.397 -3.505 -1.856 -0.703 #7.926 0.031 1.551
VI -4.57 -4.075 -4.744 -3.41 -2.462 0.329 #9.978 2.718
VII -4.152 -4.315 -4.744 -3.427 -2.901 -0.678 #4.169 #11.153
|MCSTAR critico| (alfaBonferroni=5%/49) = 3.284839
------------------------------------------
Analise de significancia pratica:
------------------------------------------
V de Cramer = 0.2404799
Grau pequeno de dependencia entre as duas variaveis nominais
Gama de Goodman-Kruskal = 0.4209053
OR generalizado = (1+0.4209053)/(1-0.4209053) = 2.453667
É praticamente o mesmo código utilizado para o exemplo de tabagismo e etilismo, com algumas diferenças:
Apenas para comparação, o teste qui-quadrado comum produziria o seguinte:
source("TCrXc_StatusOcupacional.R")
------------------------------------------
Dados:
------------------------------------------
I II III IV Va Vb VI VII
I 50 19 26 8 7 11 6 2
II 16 40 34 18 11 20 8 3
III 12 35 65 66 35 88 23 21
IV 11 20 58 110 40 183 64 32
Va 2 8 12 23 25 46 28 12
Vb 12 28 102 162 90 554 230 177
VI 0 6 19 40 21 158 143 71
VII 0 3 14 32 15 126 91 106
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: tcrxc
X-squared = 1416, df = NA, p-value = 1e-06
X^2 critico de 95% = 66.33865
Graus de liberdade (nao fornecidos por bootstrapping): 49
Heuristica de significancia (rej. H0 se X^2/gl > 2): 28.89877
------------------------------------------
Analise post hoc:
------------------------------------------
Residuos ajustados standarizados corrigidos por momento (MCSTARs):
destination
origin I II III IV Va Vb VI VII
I #28.022 #6.466 #4.851 -2.711 -0.804 -7.091 -4.336 -4.284
II #6.535 #15.194 #6.477 -0.475 0.201 -6.217 -4.43 -4.437
III 0.706 #6.01 #7.195 #3.979 2.782 -3.966 -6.129 -4.134
IV -1.369 -0.926 1.7 #6.772 0.826 0.847 -3.453 -5.132
Va -1.436 0.409 -0.87 0.701 #5.188 -1.363 0.388 -1.982
Vb -6.546 -6.397 -3.505 -1.856 -0.703 #7.926 0.031 1.551
VI -4.57 -4.075 -4.744 -3.41 -2.462 0.329 #9.978 2.718
VII -4.152 -4.315 -4.744 -3.427 -2.901 -0.678 #4.169 #11.153
|MCSTAR critico| (alfaBonferroni=5%/49) = 3.284839
------------------------------------------
Analise de significancia pratica:
------------------------------------------
V de Cramer = 0.2404799
Grau pequeno de dependencia entre as duas variaveis nominais
Gama de Goodman-Kruskal = 0.4209053
OR generalizado = (1+0.4209053)/(1-0.4209053) = 2.453667
| Especificamente para tabelas 2x2 existe o \(Q\) de Yule, dado por \[Q = {{ad-bc} \over {ad+bc}}\] Assim como o Gama de Goodman-Kruskal, varia de \(-1\) a \(+1\), mas como se aplica a tabelas 2x2, relaciona-se com o odds-ratio: \[Q = {{OR-1} \over {OR+1}}\] Uma função está implementada em psych::Yule(x, Y=FALSE), que devolve o \(Q\) de Yule quando \(x\) é uma tabela de frequências absolutas 2x2. |
Para situações em que uma variável é ordinal e a outra é nominal, por exemplo quando há um modelo de dose-resposta com um desfecho dicotômico associado a uma exposição que tem mais que dois níveis, é preferível utilizar o teste conhecido como Cochran-Armitage test for trend.
|
|
|
|
|
|
Caso ambas as variáveis sejam consideradas nominais, o teste adequado é o qui-quadrado de Pearson. Caso a variável considerada como exposição seja ordinal e o desfecho seja dicotômico (tabela \(K\)x\(2\)), o teste qui-quadrado de Cochran Armitage pode ser usado.
Para comparação, ambos os testes são executados em Teste qui-quadrado de Cochran-Armitage.R:
# Teste qui-quadrado de Cochran-Armitage.R
# Cochran-Armitage test for trend
# Menarca precoce & espessura da dobra cutânea do tríceps
# Fonte: Beckles et al. (1985) International Journal of Obesity 9:127-35,
# apud Kirkwood & Sterne (2003): Chapter 17: Chi-squared tests for 2x2 and
# larger contingency tables.
library(readxl)
library(DescTools)
# desabilita warnings
options(warn=-1)
# Tabela de Dancey&Reidy (2019, p. 275)
TCtmp <- read_excel("menarca.xlsx")
TC <- data.matrix(TCtmp)
rownames(TC) <-as.character(unlist(TCtmp[,1]))
TC <- TC[,-1]
remove(TCtmp)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
print(TC)
cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05
cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(TC,simulate.p.value=TRUE,B=1e6))
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=1), "\n", sep="")
N <- sum(res$observed)
nL <- nrow(TC) # ou dim(TC)[1]
nC <- ncol(TC) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
cat("\nTeste qui-quadrado de Cochran-Armitage:\n")
res <- DescTools::CochranArmitageTest(TC)
print(res)
x2ca <- res$statistic^2
p <- res$p.value
cat("\nX-squared trend =",x2ca, ", df = 1, ", "p-value =", p,"\n")
cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
V <- sqrt(phi2/Dim) # V ou C de Cramer
if (0 <= V & V < 0.1) {gV <- "minimo"}
if (0.1 <= V & V < 0.3) {gV <- "pequeno"}
if (0.3 <= V & V < 0.5) {gV <- "intermediario"}
if (0.5 <= V & V <= 1.0) {gV <- "grande"}
cat("\nV de Cramer =", V, "\nGrau", gV,
"de dependencia entre as duas variaveis nominais\n")
# reabilita warnings
options(warn=0)
que produz:
------------------------------------------
Dados:
------------------------------------------
menor12 maiorigual12
Pequena 15 156
Intermediaria 29 197
Grande 36 186
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: TC
X-squared = 4.7594, df = NA, p-value = 0.09404
X^2 critico de 95% = 3.841459
Graus de liberdade (nao fornecidos por bootstrapping): 2
Heuristica de significancia (rej. H0 se X^2/gl > 2): 2.379692
Teste qui-quadrado de Cochran-Armitage:
Cochran-Armitage test for trend
data: TC
Z = 2.1783, dim = 3, p-value = 0.02938
alternative hypothesis: two.sided
X-squared trend = 4.744926 , df = 1, p-value = 0.02938481
------------------------------------------
Analise de significancia pratica:
------------------------------------------
V de Cramer = 0.08768596
Grau minimo de dependencia entre as duas variaveis nominais
Compare os valores \(p\) de ambos os testes. Cochran-Armitage, ao considerar a espessura da prega cutânea como “dose”, ordinal, indica significância estatística onde o teste convencional não detectava.
Para desfecho com três ou mais categorias (os desfechos são tratados como variáveis nominais) existe função no pacote coin::chisq.test
Um exemplo está em Teste qui-quadrado de Cochran-Armitage generalizado.R:
library(coin)
library(rcompanion)
# The chisq_test function in the coin package can be used to conduct a test of
# association for a contingency table with one ordered nominal variable and
# one non-ordered nominal variable.
# The Cochran–Armitage test is a special case of this where the
# non-ordered variable has only two categories.
# The scores option is used to indicate which variable should be treated as
# ordered, and the spacing of the levels of this variable.
Job <- matrix(c(1,2,1,0, 3,3,6,1, 10,10,14,9, 6,7,12,11), 4, 4,
dimnames = list(income = c("< 15k", "15-25k", "25-40k", "> 40k"),
satisfaction = c("VeryD", "LittleD",
"ModerateS", "VeryS")))
Job <- as.table(Job)
print(coin::chisq_test(Job, scores = list("income"= c(1, 2, 3, 4))))
# Teste post hoc
print(ph <- rcompanion::pairwiseOrdinalIndependence(Job,compare="column"))
try(rcompanion::cldList(p.value ~ Comparison, data = ph, threshold = 0.05))
que produz:
Loading required package: survival
Attaching package: 'survival'
The following object is masked from 'package:epitools':
ratetable
Asymptotic Generalized Pearson Chi-Squared Test
data: satisfaction by
income (< 15k < 15-25k < 25-40k < > 40k)
chi-squared = 3.1362, df = 3, p-value = 0.3711
Comparison p.value p.adjust
1 VeryD : LittleD 0.451 0.541
2 VeryD : ModerateS 0.351 0.526
3 VeryD : VeryS 0.161 0.526
4 LittleD : ModerateS 0.698 0.698
5 LittleD : VeryS 0.242 0.526
6 ModerateS : VeryS 0.271 0.526
Este teste utiliza uma tabela tridimensional. Neste exemplo temos a contagem de mortes entre fumantes e não fumantes estratificadas por faixa etária. Esta versão de Qui-quadrado faz o teste global e, depois, estratifica pela terceira variável (faixa etária, neste exemplo).
| Faixa Etária | 18-24 | 25-34 | 35-44 | 45-54 | 55-64 | 65-74 | 75 | |
|---|---|---|---|---|---|---|---|---|
| Tabagista | Desfecho | |||||||
| Sim | Morta | 2 | 3 | 14 | 27 | 51 | 29 | 13 |
| Viva | 53 | 121 | 95 | 103 | 64 | 7 | 0 | |
| Não | Morta | 1 | 5 | 7 | 12 | 40 | 101 | 64 |
| Viva | 61 | 152 | 114 | 66 | 81 | 28 | 0 |
Os dados precisam ser reorganizados para que construamos uma tabela tridimensional requerida pela funções relacionadas a este teste. Veja como os dados foram acomodados na planila fumo_e_faixaetaria.xlsx.
Implementamos o teste em Teste Qui-quadrado de Mantel-Haenszel.R:
# Teste Qui-quadrado de Mantel-Haenszel.R
# desabilita warnings
options(warn=-1)
# Tabela de contingencia 2x2 segmentada
# Teste robusto da razao de chances (OR) de Mantel-Haenszel
# Relacao entre tabagismo e sobrevivência em 20 anos (1974-1994)
# em 1.134 mulheres adultas do Reino Unido
# Delineamento: coorte
# Fonte: APPLETON, D. R. et al. (1996) Ignoring a covariate:
# An example of Simpson's paradox. The American Statistician,
# 50(4): 340-1.
# A aspa inicial TEM que comecar na primeira coluna da linha
# e os espacamentos distintos dos dessa tabela podem causar
# problemas de na geracao da tabela horizontalizada (ftable).
library(vcd)
library(rcompanion)
library(coin)
library(readxl)
library(tidyverse)
TCtmp <- read_excel("fumo_e_faixaetaria.xlsx")
TCtmp <- gather(TCtmp, var, Contagem, -Idade)
TCtmp <- separate(TCtmp, var, c('Exposicao','Desfecho'))
TCS <- TCtmp %>% xtabs(Contagem ~ Exposicao + Desfecho + Idade, .)
print(TCS)
# Mantel-Haenszel test of the null that two nominal variables
# are conditionally independent in each stratum,
# assuming that there is no three-way interaction.
# Woolf test for homogeneity of ORs across strata on 2 x 2 x k tables:
# If significant, M-H test is not appropriate.
print(ftable(TCS)) # Display a flattened table (tabela horizontalizada)
cat("\nHomogeneidade dos ORs:\n")
print(vcd::woolf_test(TCS))
# Teste para tabela 2x2xK
cat("\nQui-quadrado de Mantel-Haenszel:\n")
print(mantelhaen.test(TCS, exact=TRUE)) # Teste exato de OR de Mantel-Haenszel
# Teste para tabela LxCxK
cat("\nTeste de Cochran-Mantel-Haenszel (coin::cmh_test):\n")
print(coin::cmh_test(TCS, distribution = approximate(nresample = 1e6)))
cat("\nAnalise estratificada por idade:\n")
print(rcompanion::groupwiseCMH(TCS,
group = 3,
fisher = TRUE,
gtest = FALSE,
chisq = FALSE,
method = "fdr",
correct = "none",
digits = 3))
cat("\n")
print(ors <- vcd::oddsratio(TCS, log=FALSE)) # Show OR for each 2x2
cat("\n")
lnors <- vcd::oddsratio(TCS, log=TRUE) # Show ln(OR) for each 2x2
print(lnors)
print(summary(lnors))
print(confint(lnors))
plot(lnors, xlab = "Faixa Etaria", ylab = "ln(OR)")
# reabilita warnings
options(warn=0)
que gera a seguinte saída:
Loading required package: grid
Attaching package: 'vcd'
The following object is masked from 'package:epitools':
oddsratio
── Attaching packages ──────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.0 ✓ purrr 0.3.3
✓ tibble 3.0.1 ✓ dplyr 0.8.5
✓ tidyr 1.0.2 ✓ stringr 1.4.0
✓ readr 1.3.1 ✓ forcats 0.5.0
── Conflicts ─────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
, , Idade = 18-24
Desfecho
Exposicao Morta Viva
NaoTabagista 1 61
Tabagista 2 53
, , Idade = 25-34
Desfecho
Exposicao Morta Viva
NaoTabagista 5 152
Tabagista 3 121
, , Idade = 35-44
Desfecho
Exposicao Morta Viva
NaoTabagista 7 114
Tabagista 14 95
, , Idade = 45-54
Desfecho
Exposicao Morta Viva
NaoTabagista 12 66
Tabagista 27 103
, , Idade = 55-64
Desfecho
Exposicao Morta Viva
NaoTabagista 40 81
Tabagista 51 64
, , Idade = 65-74
Desfecho
Exposicao Morta Viva
NaoTabagista 101 28
Tabagista 29 7
, , Idade = 75+
Desfecho
Exposicao Morta Viva
NaoTabagista 64 0
Tabagista 13 0
Idade 18-24 25-34 35-44 45-54 55-64 65-74 75+
Exposicao Desfecho
NaoTabagista Morta 1 5 7 12 40 101 64
Viva 61 152 114 66 81 28 0
Tabagista Morta 2 3 14 27 51 29 13
Viva 53 121 95 103 64 7 0
Homogeneidade dos ORs:
Woolf-test on Homogeneity of Odds Ratios (no 3-Way assoc.)
data: TCS
X-squared = 3.2061, df = 6, p-value = 0.7826
Qui-quadrado de Mantel-Haenszel:
Exact conditional test of independence in 2 x 2 x k tables
data: TCS
S = 230, p-value = 0.01591
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
0.4538410 0.9355508
sample estimates:
common odds ratio
0.6534856
Teste de Cochran-Mantel-Haenszel (coin::cmh_test):
Approximative Generalized Cochran-Mantel-Haenszel Test
data: Desfecho by
Exposicao (NaoTabagista, Tabagista)
stratified by Idade
chi-squared = 5.8443, p-value = 0.01584
Analise estratificada por idade:
Group Test p.value adj.p
1 18-24 Fisher 0.6000 1.000
2 25-34 Fisher 1.0000 1.000
3 35-44 Fisher 0.0706 0.291
4 45-54 Fisher 0.3650 0.852
5 55-64 Fisher 0.0830 0.291
6 65-74 Fisher 1.0000 1.000
7 75+ Fisher 1.0000 1.000
odds ratios for Exposicao and Desfecho by Idade
18-24 25-34 35-44 45-54 55-64 65-74 75+
0.5219512 1.2519906 0.4314109 0.7074504 0.6223718 0.9054416 4.7777778
log odds ratios for Exposicao and Desfecho by Idade
18-24 25-34 35-44 45-54 55-64 65-74
-0.65018114 0.22473479 -0.84069420 -0.34608770 -0.47421763 -0.09933253
75+
1.56397554
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
18-24 -0.650181 1.049580 -0.6195 0.53561
25-34 0.224735 0.694493 0.3236 0.74624
35-44 -0.840694 0.470642 -1.7863 0.07406 .
45-54 -0.346088 0.375584 -0.9215 0.35681
55-64 -0.474218 0.268109 -1.7687 0.07694 .
65-74 -0.099333 0.460621 -0.2156 0.82926
75+ 1.563976 2.022270 0.7734 0.43930
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2.5 % 97.5 %
18-24 -2.7073204 1.40695808
25-34 -1.1364462 1.58591573
35-44 -1.7631351 0.08174672
45-54 -1.0822181 0.39004270
55-64 -0.9997024 0.05126713
65-74 -1.0021328 0.80346776
75+ -2.3996018 5.52755287